Collective synchronization in the presence of reactive coupling and shear diversity 
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We analyze the synchronization dynamics of a model obtained from the phase reduction of the mean-field 
complex Ginzburg-Landau equation with heterogeneity. We present exact results that uncover the role of dis- 
sipative and reactive couplings on the synchronization transition when shears and natural frequencies are in- 
dependently distributed. As it occurs in the purely dissipative case, an excess of shear diversity prevents the 
onset of synchronization, but this does not hold true if coupling is purely reactive. In this case the synchroniza- 
tion threshold turns out to depend on the mean of the shear distribution, but not on all the other distribution's 
moments. 
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Reaction-diffusion systems consisting of a large number 
of degrees of freedom display a rich variety of dynamical 
regimes that are important in a wide range of fields (lHj] . In 
particular, systems composed of many interacting aggregates 
of heterogeneous, self-oscillating elements, often show oscil- 
lations at the macroscopic level as a consequence of the col- 
lective synchronization of the individual oscillators 01 0-01 • 
An appropriate model to study collective synchronization is 
the mean-field version of the complex Ginzburg-Landau equa- 
tion (CGLE) with heterogeneity, 

z.j = Zj [1 + i(uj + qj) - (1 + iqj)\zj\ 2 ] (1) 
K N 

+ +ic)^2(Zk - Zj). 

k=l 

This equation describes an ensemble of N 3> 1 globally cou- 
pled limit-cycle oscillators, each defined by a complex vari- 
able Zj = Qje l9j . Every oscillator differs from the rest in 
the natural frequency of rotation Uj and in the shear (or non- 
isochronicity) qj, which measures how the frequency of rota- 
tion depends on the oscillator's amplitude Qj. Here we con- 
sider uij and qj to be independent random variables, with a 
joint probability function p(uj, q) = g(u)h(q). 

The oscillators are coupled via a diffusive coupling of 
strength K, which has both a real (dissipative) and an imag- 
inary (reactive) components. In general a positive dissipative 
coupling drives the system to a more homogeneous state |@1 
(but see [8]). The effect of reactive coupling on synchroniza- 
tion is more intricate and strongly relies on the presence of 
shear qj ^ |@]. 

More than 30 years ago, the Kuramoto model (KM) was 
proposed as an analytically tractable system to study collec- 
tive synchronization [10]. Since then it has become a paradig- 
matic model to explain temporal organization in a large va- 
riety of natural systems far from thermodynamic equilibrium 
JUG]]]. Under some approximations, the KM can be rigor- 
ously obtained from Eq. (Q~|i. Indeed, when the mutual cou- 
pling K of the oscillators is weak, a perturbation treatment 
permits to reduce Eq. (fl]i to a set of 7Y equations for the phases 



only H, 

0j = cuj + K(qj - c) (2) 
+ KR [(1 + qjc) sin(* - 6j) - {qj - c) cos(* - 9j)\ , 

where R e 1 * = A^ 1 Yl^=i e - i6k i s tne complex order parame- 
ter. Originally, Kuramoto considered Eq. (fl} without reactive 
coupling and without shear [10]. The resulting phase equation 
(|2]i with c = qj ; = is the well-known KM. Assuming con- 
stant shear, qj = q, model (HQ) is equivalent to the so-called 
Sakaguchi-Kuramoto model Il2l IT3ll . This can be seen using 
the definition tan/3j = (qj — c)/(l + qjc), with \/3j\ < 5, 
which permits to write Eq. in the more compact form 

6 i = w * + (1 t?R )K & sin ^ - 9 i ~ + sin ft] • 

COb jjj 

As in the KM, the Sakaguchi-Kuramoto model shows a tran- 
sition from incoherence to collective synchronization at large 
enough values of K(l + qc). The synchronized solution can 
be obtained explicitly if g(uS) is a Lorentzian distribution. 

We recently reported in [ 14] that, if shear is distributed ac- 
cording to a certain probability function h(q), the onset of 
synchronization is prevented when the width of h(q) exceeds 
a precise threshold. These results were obtained assuming 
purely dissipative coupling (c = 0), and are fully analytic if 
g(oS) and h(q) are both Lorentzian. More recently IU5I1 we 
allowed ui and q to be non-independent, but still considering 
c = 0. 

Our first aim in this paper is to analyze the phase reduc- 
tion (O with reactive coupling (c 7^ 0) and independent ran- 
dom variables ojj and qj. We will show that in this case, if 
g(u) and h(q) are both Lorentzian, the onset of synchroniza- 
tion is also prevented beyond a critical value of the width of 
h(q). In the second part we address the case of purely reac- 
tive coupling, since it has physical relevance in the context 
of arrays of coupled nanomechanical oscillator s lflq - [l8ll . and 
in ion chains interacting via Coulomb forces IU9I1 . We will 
demonstrate that, in this case, the synchronization's critical 
coupling becomes fully independent of the particular shape 
of the shear distribution h(q). Finally we briefly discuss the 
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implications of this result for the KM with random coupling 
strengths recently studied by Hong and Strogatz [20, 2j]]- 

To analyze Eq. (f2]i we adopt the thermodynamic limit N —> 
oo. This allows us to define a probability density function 
(PDF) for the phases f(8, uj, q, t), such that the complex order 
parameter r = R e 1 * is 



r(t) 
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e i8 f(9,u,q,t) dOdojdq. 



The evolution of Eq. (ffjl obeys the continuity equation 

d t f= (3) 
- d e ({lj + K(q-c) + § [re-' w C(l - iq) - ex.] } /) . 

Here C = 1 + ic, and c.c. stands for complex conjugate 
of the preceding term. Next we expand / in Fourier series 
as f(0,u,q,t) = ^9{u>)h{q)YZ-oofl^^^)e m , with 
/; = f-i, an d fo = 1- Substituting the Fourier series into the 
continuity equation ([3j, we obtain the infinite set of integro- 
differential equations 



dtfi=-il[u + K(q-c)]ft 

+ M [ r * C *(l + - rC(l - iq)f l+1 ] 



(4) 



The next step is to assume that the asymptotic solutions of the 
model belong to the family of functions 



fl(u,q,t) = a{uj,q,t) 1 



(5) 



a type of ansatz discovered by Ott and Antonsen 1221 - 12411 . This 
solution of Eq. requires a to evolve according to 

d t a = — i[u) + K(q — c)]a (6) 
+ f [r*C*(l + iq) - r£7(l - iq)a 2 ] , 



with 



r*(t) 



g(ui)h(q)a(uj,q,t) dui dq. (7) 



A considerable simplification is achieved if g(oj) and h(q) 
are chosen to be Lorentzian PDFs 



<S/tt 



(w - w ) 2 + S 2 : 



7/tt 



(<? - <7o) 2 + 7 2 



(8) 



In this case the integral (0 can be solved closing the integrals 
at infinity and using the residue's theorem; notice that g(ui) — 
(27ri) _1 [(uj— Ld — iS)^ 1 — [uj— ujQ+iS)~ l ], likewise for h{q). 
The important requirement is that the complex function a can 
be analytically continued from real uj and q into the complex 
planes u = 0J r + ioJi and q = q r + iqi, inside the integration 
contours. 

It can be shown that a has no singularities in the lower half 
cj-plane |22tl. Regarding the variable q, we follow the reason- 
ing in |14|1 and find that a is analytic in the lower half complex 
g-plane for K > 0, and in the upper one for K < 0. However, 
now this holds true only if the order parameter satisfies l25ll 



R<R y 



1 



(9) 



We assume that states fulfilling this condition are correctly an- 
alyzed within this framework. As we show below, the numer- 
ical simulations fully confirm the validity of this assumption. 

Therefore, using the residue's theorem, the integrals in (0 
give 



r*(t) = a(u = oj p ,q = q p , t) = a(t), 



(10) 



where u p = ojq — iS and q p = qo T il (f° r positive and neg- 
ative K, respectively) correspond to the simple poles of the 
Lorentzian PDFs ([S]). The infinite set of ordinary differential 
equations (O then simply reduces to the single ordinary dif- 
ferential equation with complex variable 

a = -ito p a + y C(l - iq p )(l - \a\ 2 )a. 

As a = R e _i *, the equations for the order parameter inside 
the manifold defined by Eq. (0, read 

R= H + f (l + cq Tl)(l~R 2 )]R, (ID 

* = Wo + f [<? - C (lT7)](l-i? 2 ), (12) 

which we conjecture are the correct equations for the evolu- 
tion of the order parameter, as far as condition © is fulfilled. 
From Eq. (fTTT i we find that a synchronized solution bifurcates 
from incoherence at the critical coupling 



2<5 



Kr. 



which only exists if 



l+goc-7 
2(5 

l+goc+7 



if 1 
if 1 



q c > 0; 
q c < 0, 



7 < Id = |1 + qoc\. 



(13) 



(14) 



Otherwise, if 7 > 7<j, incoherence becomes the only stable 
state for all K, see Fig. [T] This result extends the one found 
in lfl4ll for c = to ensembles of oscillators globally coupled 
via both dissipative and reactive coupling. However, as it is 
depicted in Fig.[TJa), now the region of stable synchronization 
is located at positive values of K only if 



1 + q Q c> 0. 



(15) 



This inequality is the well-known Benjamin-Feir-Newell cri- 
terion for the stability of plane waves in the homogeneous 

CGLE (lUIHESIil, which is Valid in any dimension (in- 
finite in the present case). Finally, we also find that the order 

parameter of the partially synchronized solution follows 



R 2 



K-K c 
K 



with R < R x , 



(16) 



exactly as in the KM with Lorentzian g(u>) JU). However, now 
this formula holds only up to i? x , recall Eq. ©. The insets in 
Fig. [T] show numerical simulations that confirm the validity of 
Eq. ( TToT ). Remarkably, our numerical simulations indicate that 
R departs from Eq. ( [ToT l precisely when R exceeds R x . 

Unfortunately, there is no straightforward theoretical exten- 
sion of these results to more general distributions g(ui) and 
h(q), but still some reasonable conjectures can be raised: 



3 




FIG. 1. Phase diagrams with boundary d 1 3 fc for (a) qoc — —i, and 
(b) qoc — — |. Insets show the numerical results of the time averaged 
quantity R? vs. K/S, with parameters 7=4 (\K C \ = 8<5), <5 = 0.1, 
and ljo = 3. {oJj, qj}j=i,...,N were deterministically selected to 
represent the distribution ([8) with TV = 14400 oscillators. The data 
sets correspond to two different combinations of go and c, with c = 2 
(R 2 X = 1) and c = § (i?^ =4). The solid line is Eq. < TT6b . 



(i) The flip of the synchronization region when 1 + qoc 
changes sign — compare Figs.QJa) and 1(b) — is a conse- 
quence of the Benjamin-Feir-Newell criterion ( fl~5b and 
it can be supposed to be a general feature of Eq. (0. 

(ii) Although it seems difficult to prove, the divergence of 
K c at a critical value of the shear diversity [Eq. ( [Pfl il 
is likely to be a general property, as it happens in the 
purely dissipative case (c = 0) 1 14]. 

(iii) For certain distributions and parameter values, and as a 
consequence of the persistence of a fully synchronized 
solution 6j = ^ located at \K |<?(wrj) — > oo, stable syn- 
chronization and incoherence should coexist at large K 
values (as it occurs for c = and Gaussian distribu- 
tions urn ). Note that infinitesimal perturbations obey 
SBj = f (1 + <fec)[(l - N)8Qj + £^ S9i\. The Ja- 
cobian matrix has always one trivial zero eigenvalue. If 
c = and K positive, the remaining eigenvalues are 
negative, and the fixed point is stable irrespective of the 
width of h(q). However, if c ^ the fixed point be- 
comes a saddle when the qj's exceed some degree of 
heterogeneity, and hence its continuation at finite K is 
not an attractor either. In sum, under a large enough 
heterogeneity of shear, incoherence should be the only 
attractor at all K values; however, if c = 0, synchro- 
nization coexists with incoherence at large enough K — 
provided h(q) is not heavy-tailed, see Ill4ll . 



For the remainder of this article, we will concentrate on the 
case of purely reactive coupling. Motivated by the dynamics 
of nanoscale mechanical oscillator arrays, this problem was 
analyzed in detail by Cross et al. flrjl [l7tl with a coupling of 
the form ij^ J2k( z k ~ z j) m Eq. (Q~|)> and without shear diver- 
sity. To investigate the effect of shear diversity, we first write 
the phase model (|2]i without dissipative coupling. Substituting 
c = k/K in Eq. (O and letting K yields 



where k is now the total reactive coupling. Then, Eq. ( fl~3b 
suggests that in this limit the critical coupling is 



25_ 

qo ' 



(18) 



that remarkably depends on go but is independent of the 
amount of heterogeneity 7. 

The derivation of Eq. ( fl~8l ) is not rigorous because i? x = 
in this limit, and condition (0 is not fulfilled. Therefore, to 
confirm the validity of Eq. ( fl~8l ) and to determine how this re- 
sult generalizes to other distributions, next we perform the lin- 
ear stability analysis of the incoherent state of Eq. ( fT7b [ 28I1 . In 
the incoherent state all modes /;, save the trivial one /o = 1, 
vanish. The equation for the Fourier modes, related to Eq. (0]), 
is 

kI 

dtfi = -U(u - K)fi + — [r*(q - i)fi-i - r(q + i)fi+i] . 

Linearizing this equation about the incoherent state, we find 
that the only potentially unstable mode is the / = 1 



dtfi 



- «(w - k)/i 



h{J,q',t)g{u;')h{q')dui' dq' 



Let fi(ui,q,t) — b(uj,q)e xt , and neglect the trivial solution 
6 = 0. Invoking self-consistency and separating A into its real 
and imaginary parts (A = A r + iAj) yields 



2 

K 



(q — i)[X r — i(u; — k + A»)] 
A2 + (w-K + A 4 ) 2 



g(uj)h(q)duj dq. 



The interesting feature in the right hand side of this equation 
is that the integration over q is trivial and the result does not 
depend on the particular shape of h(q). Performing the limit 
A r — > + to obtain the stability threshold k c yields 



— = (qo - i) 



TTg(n c - A,) - i 



Ul — K c + Aj 



-duj 



(19) 

that only depends on h(q) through its mean value qo (defined 
as principal value if required). Finally, splitting ( fl~9l ) into its 
real and imaginary parts we obtain a system of two equations 
for the unknowns K r and A, 



(1 + qo)ng(K c - X t ) 



(l + <?o) 



_ oc to - k c + A, 



-doj 



2qo 

J 

2 



(20) 
(21) 



kR [qj s'm^ 



cos(*-0y)], (17) 



These equations can be solved for Lorentzian g(oS), and in- 
deed we recover the boundary ( fT8l ). However, note that now 
this result is stronger, since we have not imposed any con- 
straint on the shape of h(q). Here h(q) can be any distribution 
of mean qo- Additionally, an explicit value for k c can be easily 
obtained from Eq. ( f20l > if g(oj) is a uniform (top-hat) distribu- 
tion. These results for Lorentzian and uniform g(u>) are in 
agreement with those obtained in IU7I1 with h(q) = S(q — qo). 
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FIG. 2. (a) Phase diagram with boundary Q8) in the case of purely 
reactive coupling and g(ui) Lorentzian. (b) Time averages (R 2 ) ob- 
tained from numerical simulations with qo = 2, 5 = 1, N = 40000, 
and for different distributions of q: Gaussian with variance v 2 — 1, 
symmetric bidelta: h(q) = \5{q — 3) + \8{q — 1), asymmetric 
bidelta: h(q) = jS(q — 3) + \5{q + 1), Lorentzian with 7 = 1, and 
bimodal: h(q) — -, — , 1/{ \ T ,L, 2 ^ + i — , T2m2_i_i ■ m a U cases 

[l-(90-5/ 2 )] + 1 [«-(90+5/2)] ^+1 

the critical coupling is at « c /<5 = 1, as predicted by Eq. d 1 8b - 

This confirms that the phase equation is indeed a good ap- 
proximation of the amplitude equation in the limit of weak 
coupling and narrow frequency distributions. 

Figure [2 a) displays a phase diagram with the boundary 
( fT8l l. and Fig. |2fb) shows the time average (-R 2 ) obtained 
from numerical simulations using different distributions h(q) 
with common go values. As expected, in Fig. 2(b) the tran- 



sition between incoherence and synchronization occurs at the 
same value of k/5 irrespective of the distribution type. 

Finally we point out an interesting similarity between 
Eq. ([T7]l and the model recently studied by Hong and Stro- 
gatz l20l Eltl . which in our notation reads: 8j = u)j + 
qjRs'm^ — 9j). Note that here qj acts as a distributed cou- 
pling strength. Performing a stability analysis like we did 
above, we obtain that the stability border of incoherence satis- 
fies 2 = irqog(ojo), if g(u>) is unimodal and symmetric. Again, 
we obtain a formula that depends on the mean of h(q), but not 
on its shape. This result reproduces the classical Kuramoto 
relation for uniform all-to-all coupling [h(q) = S(q — qo)], 
and the critical point found in Eq. (12) of B20I1 for Lorentzian 
g(ui) and bidelta h(q). 

In conclusion, we have reported on exact results that ex- 
tend the phase models of Kuramoto and Sakaguchi ifiol [l2tl 
to situations where shear is distributed. In contrast to the re- 
cent work |3, here the coupling is not purely dissipative but 
also contains a reactive component c. In this case we also find 
that shear diversity prevents the onset of collective synchro- 
nization, but the Benjamin-Feir-Newell criterion determines 
now if the region of synchronization is located at positive or 
negative values of K. Finally, we have obtained a remark- 
able result when the coupling is purely reactive: the stability 
threshold of incoherence depends on the mean shear qo, while 
the shear diversity becomes irrelevant. 

Financial support from the Ministerio de Ciencia e Inno- 
vacion (Spain) under Projects No. FIS2009-12964-C05-05 
and No. SAF2010-16085 is acknowledged. 



[1] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence 

(Springer- Verlag, Berlin, 1984). 
[2] J. D. Murray, Mathematical biology, Interdisciplinary applied 

mathematics (Springer, New York, 2003). 
[3] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 

(1993). 

[4] A. T. Winfree, The Geometry of Biological Time (Springer, New 
York, 1980). 

[5] A. S. Pikovsky, M. G. Rosenblum, and J. Kurths, Synchroniza- 
tion, a Universal Concept in Nonlinear Sciences (Cambridge 
University Press, Cambridge, 2001). 

[6] S. H. Strogatz, Sync: The emerging science of spontaneous or- 
der. (Hyperion Press, New York, 2003). 

[7] P. C. Matthews, R. E. Mirollo, and S. H. Strogatz, Physica D 
52, 293 (1991). 

[8] H. Daido and K. Nakanishi, Phys. Rev. Lett. 96, 054101 (2006). 
[9] See 1 1, 5, 29] for the case of two coupled Stuart-Landau oscil- 
lators, and (HQIEIIII 

for the mean-field CGLE. 
[10] Y. Kuramoto, in International Symposium on Mathematical 
Problems in Theoretical Physics, Vol. 39 of Lecture Notes in 
Physics, edited by H. Araki (Springer, Berlin, 1975), pp. 420- 
422. 

[11] J. A. Acebron et al, Rev. Mod. Phys. 77, 137 (2005). 
[12] H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys. 76, 576 
(1986). 

[13] N. Nakagawa and Y. Kuramoto, Prog. Theor. Phys. 89, 313 
(1993). 



[14] E. Montbrio and D. Pazo, Phys. Rev. Lett. 106, 254101 (2011). 

[15] D. Pazo and E. Montbrio, EPL 95, 60007 (201 1). 

[16] M. C. Cross, A. Zumdieck, R. Lifshitz, and J. L. Rogers, Phys. 

Rev. Lett. 93, 224101 (2004). 
[17] M. C. Cross, J. L. Rogers, R. Lifshitz, and A. Zumdieck, Phys. 

Rev. E 73, 036205 (2006). 
[18] S.-B. Shim, M. Imboden, and P. Mohanty, Science 316, 95 

(2007). 

[19] T. E. Lee and M. C. Cross, Phys. Rev. Lett. 106, 143001 (201 1). 
[20] H. Hong and S. H. Strogatz, Phys. Rev. Lett. 106, 054102 

(2011). 

[21] H. Hong and S. H. Strogatz. larXiv: 1 108.61 171 (201 1). 
[22] E. Ott and T. M. Antonsen, Chaos 18, 0371 13 (2008). 
[23] E. Ott and T. M. Antonsen, Chaos 19, 0231 17 (2009). 
[24] E. Ott, B. R. Hunt, and T. M. Antonsen, Chaos 21, 025112 

(2011). 

[25] In the semicircular path at infinity (|g| — > oo) parametrized by 
•d,q— | q\e 1 ^ , the evolution equation of | a | evaluated at | a \ — 1 
yields at leading order dt\a\ = K(l — R\C\ cos x)\q\ sm 
where q = He -4 *, C= |C*|e ic , and X = 1>(q, t) - * (t) - C- 
Thus, the condition dt\a\ < assuring that Eq. ((5} remains 
finite, is fulfilled if R < R x = |C|~\ in ■& € (0, tt) for 
K < 0, and in ■& G (-tt, 0) for K > 0. 

[26] V. Hakim and W. J. Rappel, Phys. Rev. A 46, R7347 (1992). 

[27] I. S. Aranson and L. Kramer, Rev. Mod. Phys. 74, 99 (2002). 

[28] S. H. Strogatz and R. E. Mirollo, J. Stat. Phys. 63, 613 (1991). 

[29] D. G. Aronson, G. B. Ermentrout, and N. Kopell, Physica D 41, 



403 (1990). 



